## Simulations of the BSC and Rep(3)

from PyM import rd_float, hd, show
from PyM import factorial

# This function simulates a binary symmetric channel 
# of bit-error probability p when a binary message x 
# is sent through it. Here x is assumed to be a 
# streem of bits and 0<=p<=1
def BSC(x,p):
    y = ''
    for b in x:
        r = rd_float()
        if r<=p: 
            if b=='0': y += '1'
            else: y += '0'
        else: y += b
    return y

# Encoding and decoding Rep(3)
def rep3enc(u):
    x = ''
    for b in u:
        x += 3*b
    return x

def rep3dec(y):
    u = ''
    while len(y)>2:
        c = y[:3]; y = y[3:]
        if c.count('1')>1:
            u += '1'
        else: u += '0'
    if y=='':
        return u
    print('rep3dec stopped at undecodable tail', y)
    return u

# Testing BSC together with Rep(3) coding and decoding
def rep3test(u,p):
    x = rep3enc(u)
    y = BSC(x,p)
    v = rep3dec(y)
    E = len(u)*p   # expected number of errors with no coding
    Ec = len(u)*(3*p**2*(1-p)+p**3)  # " " " " with coding 
    m = hd(u,v)                     # number of errors with coding
    return (E, round(Ec,3) ,m)

# Example
L = 100000; p = 0.005
# E = L*p = 1500, EE = 22.425
# The value of N, the number of errors with coding, 
# will fluctuate about EE when the test is repeated.  
# For higher values of L, say 500000 or higher, 
# the computation may take too much time.
show(rep3test(L*'0',p))

## Remarks
'''
Under the conditions of the above example, the probability p_m of getting m errors in a trial is 
     binom{L}{m}{p'}^m(1-p')^{L-m}, p'=3p^2(1-p)+p^3. 
     
This 'Bernoulli expression' can be well approximated by the 'Poisson distribution' 

     exp{-lambda}  lambda^m/m!, \lambda = L*p'= 7.475, 

which means that the p_m are proportional to lambda^m/m! with scaling factor
K=e^{lambda}.

Next computation shows that the p_m increase from $m=0$ up to $m=7$ and then decrease. Moreover, the chances that m <= 7 are slightly higher than the chances of getting m>7
(the odds are 53 to 47). The figures also show that the values of m are rather dispersed
around m = 7. For instance, the value of m falls 

     42% of times within the interval $[6,8]$,
     65% within $[5,9]$,
     80% within $[4,10]$, 
     90% within $[3,11], and so on. 
'''

l = 7.475
P = [(m,l**m/factorial(m)) for m in range(19)]
P = [(p[0],round(p[1],1)) for p in P]
print(P)

from math import exp
K = round(exp(l),1)
print(K)

K_low = round(sum([P[m][1] for m in range(8)]),1)

print(K_low)

p_low = round(K_low/K,2)

print(p_low)

print(1-p_low)

#
dm = 2
K_around= round(sum([P[m][1] for m in range(7-dm,8+dm)]),1)
print(K_around)
P_around = round(K_around/K,2)
print(P_around)